function [R,F,Ab,LK] = mesh_center_2d_resp(keff, in)
% TWO-DIMENSIONAL, MESH-CENTERED, MULTIGROUP DIFFUSION
% response matrix generator

dat   = in.dat;
%abdat = in.abdat;
numg  = in.numg;
numm  = in.numm;
%numel = in.numel;
order = in.order;

R    = zeros( in.numtypes, 4*numg*(order+1), 4*numg*(order+1) );
Rinc = zeros( numg*(order+1), numg*(order+1) );
Rfar = Rinc; Rlef=Rinc; Rrig=Rinc;
F    = zeros( in.numtypes, 4*numg*(order+1) );
Ab   = F;
tmpF = zeros( 1, numg*(order+1) );
tmpAb = tmpF;
LK = zeros( in.numtypes, 4,  numg*(order+1) );


% loop over all element types ---------------------------------------------
for el = 1:in.numtypes
    
    mt  = in.mt(:,:,el);
    
    BCL = 0;
    BCR = 0;
    BCB = 2;
    BCT = 0;
    
    src(:,:,1) = [ 0
        0 ];
    src(:,:,2) = [ 0
        0 ];
    
    % loop over legendre order and group ----------------------------------
    for RG = 1:numg
        for RO = 0:order
            
            %---- GENERATE THE COEFFICIENTS -------------------------------
            [AB, AL, AC, AR, AT, Si, Sh, siga, dx, dy, dv, dc, sct, chi, alphaB] = ...
                mesh_center_2d_coef_resp( dat, numg, numm, in.xcm, ...
                in.xfm, in.ycm, in.yfm, mt, src, 1,[BCL BCR BCB BCT], RO, RG );
            
            %---- SOLVE----------------------------------------------------
            N = sum(in.xfm); M = sum(in.yfm);
            % pad end of AB,AL
            AB = [AB' zeros(numg,N)]';
            AL = [AL' zeros(numg,1)]';
            % pad top of AR,AT
            AT = [zeros(numg,N) AT']';
            AR = [zeros(numg,1) AR']';
            phi = zeros( N*M, numg );
            
            % build cell structure of sparse matrices
            for g = 1:numg
                A{g} = spdiags([AB(:,g) AL(:,g) AC(:,g) AR(:,g) AT(:,g)], ...
                    [-N -1 0 1 N], N*M,N*M);
            end
            
            nusig = Sh;  % just rename it right away to avoid confusion
            s = 0*ones(N*M,1)/(in.ycm(end)+in.xcm(end));  % uniform guess
            phieps = 1e-9;  mxit = 2000;  i= 0; phierr = 1;
            while (phierr > phieps && i < mxit)
                phi0 = phi;
                phi(:,1) = A{1} \ ( Si(:,1) + chi(:,1).*s(:,1)/keff ) ;
                for g = 2:numg
                    scsrc = zeros(N*M,1); % reset
                    for gg = 1:numg
                        scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg); %<-----
                    end
                    phi(:,g) = A{g} \ (chi(:,g).*s(:,1)/keff + scsrc(:,1) + Si(:,g)); 
                end
                phierr = norm( phi0 - phi );
                s(:,1) = 0; % reset fission source
                for g = 1:numg
                    s(:,1) = s(:,1) + Sh(:,g).*phi(:,g);
                end 
                i = i + 1;
            end
            if (i==mxit)
                disp('exceeded mxit...')
                disp(phierr)
                disp('--')
            end
            
            %==============================================================
            
            % reponses due to  RG,  RO incident current
            [CRL CRR CRB CRT LL LR LB LT] = ...
                cresp_mc(phi,dx,dy,dc,mt,in.xfm,in.yfm,numg,1*alphaB/2);
            tmp1=0; tmp2=0;
            for g=1:numg
                tmp1 = tmp1 + sum( phi(:,g).*nusig(:,g) );
                tmp2 = tmp2 + sum( phi(:,g).*siga(:,g) );
            end
            tmpF(1, (RO)*numg+RG) = tmp1;
            tmpAb(1, (RO)*numg+RG) = tmp2;
            % compute leakage
            LK(el,1,(RO)*numg+RG) =  LB; %incident
            LK(el,2,(RO)*numg+RG) =  LT; %far
            LK(el,3,(RO)*numg+RG) =  LR; %to right of incident
            LK(el,4,(RO)*numg+RG) =  LL; %to left of incident
            
            % determine legendre coefficients -----------------------------
            % Lp's should be put in a function and moved outside this loop
            % determine legendre polynomial for incoming current
            NN = sum(in.xfm)-1;    % max degree
            k = 0:NN;         % interval
            maxord = order;  % maxorder to which we compute P's
            P = zeros(NN+1,NN+1);
            P(1,1:NN+1) = 1;
            P(2,1:NN+1) = 1 - (2*k)/NN;
            for j = 3:NN+1 %otherwise N-1
                i = j-2;
                P(i+2,1:NN+1) = ( (2*i+1)*(NN-2*k).*P(i+1,:) - i*(NN+i+1)*P(i,:) ) ./ ...
                    ((i+1)*(NN-i)) ;
            end
            for i = 0:maxord
                a = (2*i+1);
                b = factorial(NN+i+1)/factorial(NN);
                c = factorial(NN)/factorial(NN-i);
                w(i+1)=a*c/b;
            end
            %x = xint;
            
            for G=1:numg
                % f def'n ASSUMES four-way symmetry
                f(1,:) = CRB(end:-1:1,G)';
                f(2,:) = CRT(:,G)';
                f(3,:) = CRR(end:-1:1,G)';
                f(4,:) = CRL(:,G)';
                
                % now we want to expand f into fi
                for ord=1:maxord+1
                    ff(:,ord) = w(ord) * f(:,:) * P(ord,:)';
                end
                %side=1;
                Rinc( G:numg:end, (RO)*numg+RG ) = ff(1,:)';
                Rfar( G:numg:end, (RO)*numg+RG ) = ff(2,:)';
                Rrig( G:numg:end, (RO)*numg+RG ) = ff(3,:)';
                Rlef( G:numg:end, (RO)*numg+RG ) = ff(4,:)'; 
            end
            
            % end loop over legendre order and group ----------------------
        end
    end
    % ---------------------------------------------------------------------
    R(el,:,:) = [  Rinc  Rfar  Rrig  Rlef
                   Rfar  Rinc  Rlef  Rrig
                   Rlef  Rrig  Rinc  Rfar
                   Rrig  Rlef  Rfar  Rinc ];
    F(el,:)   = [ tmpF tmpF tmpF tmpF];
    Ab(el,:)  = [ tmpAb tmpAb tmpAb tmpAb];
    
    % end loop over element types -----------------------------------------
    
    
end
